function [w,SIG] = l2relax0(y, X, tau, opt)
% L2-relaxation estimation (by the default "fmincon" package). 
% opt = 1. no shrinkage
%       2. linear shrinkage
%       3. nonlinear shrinkage

    % 0. house clean
    N = size(X,2);
    if isempty(y)        
        E = X;              % if no input y
    else
        E = repmat(y,1,N) - X;            
    end
    E = E-repmat(mean(E,1),size(E,1),1);            % demean error VC
    % 1. setup parameter
    SIG = SIG_fun(E,opt);
    a0 = zeros(N,1);
    % 2. compute for a
    options = optimset('display','none');
    a = fmincon(@(a) opt_fun(a,SIG,tau),a0,[],[],ones(1,N),0,[],[],[],options);
    % 3. estimate w
    A = (eye(N)-ones(N)/N)*SIG;
    w = A*a + ones(N,1)/N;
end

function R = opt_fun(a,SIG,tau)
    N = length(a);
    A = (eye(N)-ones(N)/N)*SIG;
    R = a'*(A'*A)*a/2 + ones(1,N)*SIG*a/N + tau*sum(abs(a));
end

function SIG = SIG_fun(E,opt)
    if opt == 1         % no shrinkage
        SIG = cov(E);
    elseif opt == 2     % linear shrinkage following Ledoit and Wolf (2004) 
        SIG = cov1para(E);
    elseif opt == 3     % nonlinear shrinkage following Ledoit and Wolf (2020)         
        SIG = analytical_shrinkage(E);
    end
    if isnan(sum(sum(SIG)))
        SIG = cov(E);        
    elseif isinf(sum(sum(SIG)))
        SIG = cov(E);        
    end    
end




